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Abstract 

We report a simulation study of the gas-liquid critical point for the square- well potential, for val- 
ues of well width 6 as small as 0.005 times the particle diameter a. For small 5, the reduced second 
virial coefficient at the critical point Sl'^ is found to depend linearly on 6. The observed weak linear 
dependence is not sufficient to produce any significant observable effect if the critical temperature 
Tc is estimated via a constant -B^'^ assumption, due to the highly non linear transformation between 
-62'^ and Tc- This explains the previously observed validity of the law of corresponding states. The 
critical density pc is also found to be constant when measured in units of the cubed average distance 
between two bonded particles (1 + 0.55)cr. The possibility of describing the 5^0 dependence with 
precise functional forms provides improved acccurate estimates of the critical parameters of the 
adhesive hard-sphere AHS model. 

PACS numbers: 05.20.Jj, 61.20.Ja, 64.70.Ja 
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I. INTRODUCTION 



Investigation of protein and colloidal systems has focused the attention of the scientific 
community on the phase-diagram behavior of short-range attractive potentials and of the 
role of the range of interaction in controlling the thermodynamic and dynamic properties of 
the system Q, 3,30,3. Colloidal particles, due to their nano or microscopic size are often 
characterized by effective interactions 61] whose range is significantly smaller than the particle 
diameter. Under these conditions, it has been argued that the actual shape of the potential is 

n n 

irrelevant and that the thermodynamics[7|], as well as the dynamics [8|], of different systems 
approximately satisfy an extended .aw :;„d,„gstate.B. This law allows for a 
comparison between different systems, once the effective diameter of the particle is known 
(i.e., when the repulsive part of the interaction can be mapped into an equivalent hard- 
sphere diameter 1Q|]). It has been proposed that the second virial coefficient, normalized 
by the corresponding hard-sphere second virial coefficient, B2 may act as a proper scaling 
variable. Therefore systems with equal second virial coefficient and effective diameter should 
have similar thermodynamical properties. 

The adhesive hard-sphere (AHS) potential ll| , the limiting behavior of an infinitesimal 
interaction range coupled to infinite interaction strength such that B2 is finite, has also 
received significant attention. For this potential B2 = I — l/4r, where r, which acts as an 
effective scaled temperature, is the so-called stickiness parameter. Despite the known ther- 
modynamic anomalies 12] , an analytic evaluation of the (metastable) critical point within 
the Percus-Yevick closure with both the energy and the compressibility routes for this po- 
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tential is available [ill. . In the energy route the critical point is located at B2 
(tc = 0.1185) and number density p = 0.609. 

The availability of analytic predictions for this model has favored its application in the 
interpretation of experimental data for several disparate colloid (and protein) systems 3, ^, 
Q, S, 14, [15I, [16 1, an apphcation whose validity has been reinforced by the extended law 
of corresponding states. For this reason, it is important to try to accurately estimate the 
properties of the AHS model as a reference, to support existing predictions or to suggest 
improvements to available theoretical approaches. Numerical simulations of the AHS model 



have been attempted in the past 



13, 



18 



19|. A recent effort in the direction of evaluating 
the phase diagram of the model has been provided by Miller and Frenkelj2o|. based on an 
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ingenious identification of tlie appropriate Monte Carlo (MC) moves for this potential 17l.ll8|. 
Their study provides an estimate of the location of the critical point at = —1.21(1) 
(re = 0.1133(5)) and pc = 0.508(10). 

In this article we propose a different approach to the numerical evaluation of the criti- 
cal properties of the AHS model, based on extrapolation of standard grand canonical MC 
simulation results for a sequence of square well (SW) potentials with progressively smaller 
attraction ranges 6 (down to 6 = 0.005, in units of the hard-sphere diameter a). 

The SW potential is defined as: 

oo if r < o" 
— e if cr < r < (T + (5(7 
if r > (J -|- (5(7 

The SW fluid has been profusely studied 
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been shown 
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27| for (5 > 0.1. It has 



25( 1 that for b < 0.25 gas-liquid separation becomes metastable with respect 



to the fluid-solid equilibrium. Despite its metastable character, investigation of smaller 8 
values retains its importance, since the crystallization time is often much longer than the 
experimental one and gas-liquid phase separation is readily accessed (an effect facilitated by 
the large difference between the fluid density and the crystal density and, in experiments, 
by the intrinsic sample polydispersity) . 

Despite the importance of the SW model in relation to the AHS potential, no studies 
of the (5-dependence of the critical point location has been previously reported for very 
small (5. This is in large part due to the fact that for smaller and smaller 5, the critical 
temperature signiflcantly decreases. Indeed, according to the constant prediction of 
Noro and Frenkelj?!, it should vary as 



;i + 5)3-i;j 

making bond-breaking (changes of the particle energy of order e) events rarer and rarer in 
the simulation. Moreover, the size of the translational step in the MC code is of the order of 
(5. On the other hand, the location of the critical point becomes more and more metastable 
which, in principle, poses a limit to the smallest b which can be studied. 

Despite these numerical difficulties, we have been able to estimate the location of the 
critical point down to (5 = 0.005(T. We present here the b dependence of the critical temper- 
ature and density and the values of the second virial coefficient and energy at the critical 



point. In all short linear extrapolation to 5 = provides novel accurate estimates 

of the corresponding quantities for the AHS model. 



II. MONTE CARLO SIMULATIONS 

We have simulated the SW system in the grand canonical (GC) ensemble in order to 
locate the gas-liquid critical point for different 6 values. The critical point is identified by 
mapping the grand canonical density distribution onto the universal Ising model distribution, 



281]. Histogram re wight ing [29] was 



following the method described by Bruce and Wilding 
used to achieve an accurate estimate of the critical point, and the field mixing parameter 
was always found to be negligible. We define a Monte Carlo step as one hundred trial moves, 
with an average of 95% translation and 5% trial insertions or removals of one particle in 
the system. A translational move is defined as a displacement in a random direction by 
a random amount between ±5/2. We simulate different realizations of the same system 
(at the same chemical potential and temperature) to improve statistics. Simulations lasted 
more than 10^ MC steps. We have studied cubic boxes of side 5a and/or 8a, to estimate the 
importance of finite-size effects. Additionally, for 6 = 0.05 we have studied several other box 
sizes to estimate the deviation of the value of the critical parameter for the bulk limit case. 
Proper finite size studies for 6 < 0.05 are at this moment numerically prohibitive. Histogram 



reweighting was used to map the density distribution of the fluid onto the universal critical 



order parameter distribution of the Ising model [29|], thereby reaching an accurate estimate 
of the critical point. The field mixing parameter was always been found to be negligible. 

We have occasionally observed a transition to a more dense stable phase, signaling that 
indeed the values of the chemical potential studied admit metastable fluid solutions. In 
all cases where a transition to a more dense phase was observed, the simulation was inter- 
rupted and the 10% of configurations saved just before the transition were disregarded. This 
procedure is shown graphically in Fig.-[T1 
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FIG. 1: Density evolution for a single run of a SW fluid with 5=0.01, T* = 0.233, and ^/kT = 
—2.41 (close to the critical point). At long time, a transition to a dense phase is observed and, 
consequently the simulation on the right side of the vertical line is disregarded. 

III. RESULTS 

Fig.-l2]-(a) and (b) show the 5 dependence of the critical temperature and the correspond- 
ing critical second virial coefficient B2^. For SW, the virial can be calculated as 

B;{T) = l-{{5 + lf-l){e^'^ -I) (3) 

For small 5, a clear linear dependence of i?2'^ is observed. The data are well represented 
by a functional form B2^{5) = —1.174 — 1.774 6. The extrapolated value of -82'^ for 6 = 
(-82^ = —1.174) is slightly higher than the value —1.21(1) estimated by Miller and Frenkel for 
the AHS potential. The 6 dependence of -82^^ formally violates the idea that i^g^^ is the correct 
scaling variable for collapsing the phase-diagram of different short-range attractive potentials 
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onto a single master curve. Nevertheless, the constant ^2*^ approximation is sufficiently good 
to explain the Tc dependence, due to the non-linearity of the transformation (Eq. [2]). To 
prove this point we show in Fig.-[21-(a) predicted according to -82'^ = —1.174. In this 
representation, the assumption of constant i?2^ is sufficient to describe the 6 dependence of 
Tc up to 5 = 0.05, with an error less than 1% (growing with 6). 

Next we compare the energy per particle of the system at the critical point in Fig.-[2]-(c), 
to provide a measure of the number of contacts per particle. This quantity also shows a 
linear dependence on S. 

The study of the 6 dependence of the critical density pc has received considerably less 
attention than the 6 dependence of T^. Recently Ref. [30| suggested a plausible relation for 
the 6 dependence of pc in the limit of small well width: 

The relation is based on the hypothesis that pc should be constant if measured using the 
average distance between two bonded particles {l + 6/2)a as the unit of length. The relation 
was also supported by a potential energy landscape interpretation of the generalized law of 



corresponding states 



311 ] which shows that configurations with the same Boltzmann weight 



are generated under an isotropic scaling (to change the inter-particle distances preserving 
the same bonding pattern) and a simultaneous change of both 6 and T such that the bond 
free-energy remains constant. 

Figure [3]-(a) shows the calculated evolution of the critical density with the range of the 
interaction. It also reports previous estimates for the same system 4[] as well as the critical 
density for the AHS model from Ref. [2^ . As the range of the SW potential is reduced, the 
critical density becomes higher, since a higher local density is required to generate bonded 
configurations. Data for 6 < 0.1a are properly represented by Eq. HI with a resulting fitting 
parameter p(0) = 0.552. This value is significantly higher than the AHS critical density 
Pc = 0.508 reported in Ref.-^. 

For completeness we show in Figure [3]-(b) the 6 dependence of pc/ksTc, where pc is the 
value of the chemical potential at the critical point. pc/ksTc also shows a linear dependence, 
with an intercept at —2.394, corresponding to a critical activity exp{pc/kBTc) = 0.091, to 
be compared with the corresponding value of 0.087 of Miller and Frenkel. 

It is known that the finite size of the system modifies the position of the critical point 32] . 
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FIG. 2: Dependence of the critical temperature (a), of the second virial coefficient (b) and of the 
potential energy (c) on the range of the potential 6. Circles and crosses label simulation data of 
this work with box sides 5cr and 8a respectively. The line in (a) is the theoretical prediction for 
the critical temperature provided by EqlJl The line in (b) corresponds to the best linear fit of the 
simulation data for 5 < 0.10 (—1.174 — 1.7745). Open squares correspond to the simulation data 



from Ref. 



J]. Filled squares are the AHS i?!'^ result from Ref. 



201]. The inset presents the whole 



range of S values studied. 




FIG. 3: Evolution of the critical density (a), and chemical potential (b) with the range of the 
potential 6. Circles and crosses correspond to box sides L = 5a and L = Sa respectively. Empty 
squares in (a) are data from Ref. [4]. The lines correspond to the best fit of (a) pc for 6 < 0.10 
(0.5519/(1 +(5/2)^^and (b) Hc/kBTc for 5 < 0.10 (-2.394 - 1.431(5). Filled squares represent AHS 
results from Ref. [20,]. The inset presents the whole range of 6 values studied in this work. 
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A precise estimate of the critical point requires a complete finite-size scaling study to ex- 
trapolate the results to an infinite system. We have not attempted to perform such a careful 
study since it would be computationally prohibitive for the small ranges studied here and 
we have limited ourselves to four different box sizes (L = 5, 6, 8 and 10) for 6 = 0.05 only. 
The results are reported in table [11 As already suggested by the minor differences in the 
L = 5 and L = 8 data shown in the previous figures, no significant changes in the critical 
parameters are observed. From the scatter in the data (similar to that found by Miller and 



Frenkelj20|) it is possible to estimate errors in the critical parameters. The resulting values 
of and errors in the critical parameters are Tc = 0.3658 ± 0.0005, pc = 0.513 ± 0.008 and 
fic = 0.9064 ± 0.0008. These values of the critical temperature and density suggest that the 
Percus-Yevick energy rout e I13I is even better than previously thought, despite the fact that 
the compressibility route ll|] results seem to be more often used and cited. 



TABLE I: Critical parameters for the SW system of 6 = 0.05 obtained with four different boxes of 
side L = 5, 6, 8, 10. 



L 


Tc 


pc 


Pc 


5 


0.3660 


0.516 


-0.9042 


6 


0.3660 


0.516 


-0.9051 


8 


0.3658 


0.513 


-0.9062 


10 


0.3657 


0.511 


-0.9069 



IV. DISCUSSION 

Our study provides a set of values for the limiting AHS case, based on an accurate 
extrapolation of the critical parameters of the SW potential to 5 — » 0. These values are 
outside the error bars of Miller and Frenkel's investigation. In particular, both the critical 
density and the critical virial appear to be higher than the previous estimates. 



The special techniques employed in the simulation of the AHS system [18|, l20|, 1331] only 
consider moves that make or break up to three contacts. A particle can readily gain more 
than three contacts, since higher coordination states are established by a succession of such 
moves. Apparently, however, the constraint on the possible moves disfavors the formation 



of small nuclei of solid phases, since crystallization was extremely rarely observed with this 
algorithm, and then only in fluids with a very high mean reduced density (greater than 0.9). 
In the SW simulation, the transient solid-like nuclei are more readily formed (and indeed 
we do occasionally observe crystallization during the simulation) suggesting that the SW 
simulations sample a larger region of configuration space than that accessible with the AHS 
algorithm. This could indeed explain why the critical density extrapolated from the SW 
simulations is significantly higher than that calculated previously. 

To support this interpretation we have compared the distribution of the number of con- 
tacts per particle (proportional to the energy of the particle) for the AHS and a SW with 
6 = 0.01 at the same virial coefficient (slightly above the critical one) and same density 
for three different state points. The results of MC simulations in the NVT canonical en- 
semble, for different densities are reported in Fig. |H The distributions of the number of 
contacts per particle are coincident for low densities, but discrepancies appear as the den- 
sity is increased, confirming that the algorithm used in Ref. 201] explores configurations 
with a somewhat smaller coordination number than does the standard MC SW simulation. 
Figure H] nevertheless confirms that the AHS algorithm permits the formation of high co- 
ordination states. To avoid any potential artifact due to the a priori unknown mapping in 
the density between the AHS and the SW potential, we have also repeated the calculation 
at a lower density, scaled according to Eq. HI However, as shown in Fig. |U even when the 
density is scaled to account for the different bond distance in the AHS and SW models, at 
high density the disagreement between the two set of simulations remains. 



V. CONCLUSIONS 



We have reported a simulation study aimed at evaluating the dependence of the critical 
parameters of the short-range square well potential for interaction ranges approaching zero, 
down to 6 = 0.005(7. The resulting values for 52*^ and /3c/ic in the range 0.005 < 6 < 0.1 are 
very well described by a linear dependence on 6, providing an estimate for the 5 — > limit. 
From the resulting value of i?2^ at 5 = 0, it is possible to evaluate also the corresponding 
critical value of the AHS model via the relation B2 = 1 — l/4r. In the same range, the 
critical density is well represented by the previously proposed functional form of Eq. |H More 
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FIG. 4: Probability distribution of the number of contacts per particle for the AHS system and the 
SW with 6 = 0.01, for different densities. The comparison has been made at a slightly supercritical 
stickiness parameter r=0.120, in the case of the AHS, and the corresponding temperature provided 
by Eq. [2]for the SW of 5 = 0.01. The simulations were performed in the NVT canonical ensemble 
with volume V = 512a^. Symbols correspond to AHS simulation data, for />=0.195 (circles), 0.488 
(squares), 0.781 (diamonds). Dashed lines correspond to SW simulation data and dotted lines to 
SW simulations at the rescaled densities (see Eq. H]) 0.192, 0.481, 0.770, respectively. 

precisely, our extrapolation for the AHS limit is: 



PcficiS = 0) 



BTiS = 0) 



p^{S = 0) = 0.552 ±0.001 



Tc = 0.1150 ±0.0001 



-2.394 ±0.001 



1.174 ±0.002 



(5) 
(6) 
(7) 
(8) 
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TABLE II: Critical point parameters for the SW fluid for width 5, resulting from simulations of 
boxes of side L = 5 (top part of the table) and 8 a (bottom part of the table) respectively. 



6 


^ c 


Pc 


Pc 


0.005 


0.2007 


0.542 


-0.4817 


0.01 


0.2328 


0.540 


-0.5614 


0.02 


0.2769 


0.538 


-0.6693 


0.03 


0.3106 


0.530 


-0.7575 


0.04 


0.3398 


0.522 


-0.8333 


0.05 


0.3660 


0.516 


-0.9042 


0.10 


0.4780 


0.478 


-1.2120 


0.05 


0.3658 


0.513 


-0.9062 


0.10 


0.4780 


0.478 


-1.2138 


0.20 


0.667 


0.421 


-1.7812 


0.30 


0.847 


0.376 


-2.3505 


0.40 


1.029 


0.339 


-2.9525 


0.50 


1.220 


0.310 


-3.6002 


0.60 


1.430 


0.287 


-4.3051 


0.70 


1.665 


0.272 


-5.0890 


0.80 


1.940 


0.263 


-5.9636 



where the error bars are based on the spread of results from the different box sizes in this 
study. 

Three important observations are in order: 

i) Even for short-range interactions, -Bg*^ shows a shght linear dependence on 5, apparently 
contradicting the law of corresponding states. While this is technically correct, one must 
also remember that the small changes of -Bg'^ with 6, in the range < 5 < 0.05 are not 
sufficient to produce any significant observable effect in the critical temperature estimated 
using a constant -Bg'^ assumption via Eq. O This is due to the highly non-linear relationship 
between the two quantities, explaining the success of the law of corresponding states in 
the interpretation of simulation and experimental data of short-ranged attractive potentials. 
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From a more academic point of view, the law of corresponding states may still hold but 
with a scaling variable more complicated than the reduced virial coefficient itself and with 
a proper scaling of the density. 

ii) The "best" critical parameters of the 6 = case are found to be different from those 
reported by Miller and Frenkel |2Ci] . We believe this discrepancy is related to an incomplete 
mapping of the configuration space which manifests itself in a less complete sampling of the 
dense region. The extreme rarity of any hint of crystallization in the simulations of Miller 
and Frenkel is consistent with this proposed explanation. 

iii) The higher critical temperature and density established by the SW extrapolation 
suggest that the Percus-Yevick energy route offers a more accurate estimate of the AHS 
critical point than the Percus-Yevick compressibility route. 

We acknowledge support from MCRTN-CT-2003-504712 and J. L. MERG-CT-2006- 
046453. We thank D. Frenkel and G. Foffi for helpful discussions. 
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